09 - Remote Sensing with Landsat package
Working with satellite imagery in R
In this exercise you learn to use the functions of the
landsat package to process satellite imagery, specifically,
how to:
- streamline histograms so you can combine images
- manipulate multi-band imagery,
- extract and create new data from images, such as NDVI and SAVI
- classify image values using
kmeansunsupervised algorithm to detect similar areas - segment features in satellite imagery
Task 1: Set up your workspace
Start by installing the required packages first.
rasterVis and RColoBrewer are used for
visualisation of rasters, lattice and
latticeExtra for extra graphical utilities, while
landsat provides the imagery samples and rgl
allows for 3D visualisation.
Task 2: Pre-processing of Landsat datasets
Landsat packages offers sample Landsat satellite imagery decomposed
into single bands, and labelled by the month of capture (25 November
2002). You will load and practice raster analysis with these 300x300
pixel samples. Note that each single band image shows as panchromatic
(greyscale) with values ranging from 0 (white) to 255 (black), with the
interim colors being shades of grey. This is unsigned 8 bit imagery. The
number behind the filename (e.g. nov3) refers to the color
band of the image: 1 = Green, 2 = Blue, 3= Red, 4 = Near-infrared. Once
you combine these color bands and plot with plotRGB(), you
can see the true- or false-color imagery depending on which color you
load into which RGB channel.
The images are in the SpatialGridDataFrame format and
need to be converted to raster format before
manipulation.
library(landsat)
## load indvidual band image data from landsat package
`?`(nov)
# load band#3 red channel of the image
data(nov3)
plot(nov3)Task 3: Load elevation data and plot it in 3D
dem in the landsat package denotes a ‘digital elevation
model’. Once you load it you can convert it to Formal Raster Layer with
raster() function and then continue processing using the
usual raster package functions.
plot3D() is a neat function in the
rasterVis package, which opens a separate windows and plots
elevation values in 3D space if present in the raster. You need to close
the window if you want to update or plot another object.
Plot in 3D
Task 4: Load and explore RGB data components
# let's load data for July and explore it
data("july1")
data("july2")
data("july3")
data("july4")
j1 <- raster(july1) # blue
j2 <- raster(july2) # green
j3 <- raster(july3) # red
j4 <- raster(july4) # near-infrared
## check out the image histogram\t
plot(j1)Task 5: Plot RGB image
# take the June data and create an RGB image and drape it over the 3D model
### Reorder to R - G - B and create a multi-layer rasterBrick and also a
### false-color rasterStack!
myRGB <- brick(j3, j2, j1) # brick creates new object
myCIR <- stack(j4, j3, j2) # stack stores connections only
### let's see how the NIR, R, and G bands relate (from lattice)
splom(myCIR, varname.cex = 1) # scatter plot matrix
Ok, you can probably see something?
Task 6: Manipulate image rendering by histogram stretching
First, let’s use histogram stretch
Next, let’s try linear stretch
Any idea what the red color represents in myCIR?
Finally, drape one of the images over a 3D model. This bit can be a
bit touchy and take a while to get to work. I tend to run the
plot3D() lines alone to generate the 3D view in a pop-up
window rather than rmarkdown output. Everytime you wish to refresh the
view, you must close the pop-up window.
## finally...
rglwidget() # this widget helps get the first view rendered in rmarkdown, refreshing however is more tedious
`?`(rglwidget())
t <- plot3D(dem, col = rainbow(255)) ## you need to close RGL device manually first and then run this line!
# Save it to a file. This requires pandoc
filename <- tempfile(fileext = ".html")
htmlwidgets::saveWidget(rglwidget(), filename)
browseURL(filename)
plot3D(dem, drape = j4) ## should drape image j4 over DEM, if problematic try in .R script and watch for a pop-up widget. Ask Adela to demo!Additionally, the histograms might need some adjustment to balance/equalize them before plotting. Notice how the November and July value scales are different at 80 and 250 max respectively. Clearly there is much more variation in July growth and we need to stretch the November histogram a bit to have a comparable view.
# histogram/color adjustments
data(nov3)
data(july3)
par(mfrow = c(2, 2))
image(nov3, main = "nov")
image(july3, main = "july")
plot(nov3, main = "nov")# let's create a matching histogram by stretching nov3 to july3
nov3.newH <- histmatch(master = july3, tofix = nov3)
image(nov3.newH$newimage, main = "new nov") # look at it on the fly
# write the new image with a new histogram
n3new <- raster(nov3.newH$newimage)
# convert existing images to raster format before plotting
n3 <- raster(nov3)
j3 <- raster(july3)
# plot (the mfrow is set up to work in R script, so just click through the
# rmarkdown view)
par(mfrow = c(1, 3))#### most important corrections are atmospheric and topographic however, these
#### are too complex to cover here...see package help
par(mfrow = c(1, 1)) # remember to change plotting to a single window if in .R script
# Plot the equalized-histogram images
plot(raster(july3), main = "july")Well done on equalizing histograms over the same area. Now you can see different kinds of vegetation thriving at different times of year.
Task 7: Create new information from satellite imagery: NDVI
Let us calculate the Normalized Difference Vegetation Index (NDVI) and see where the vegetation grows most in our Landsat image: Remember the formula for NDVI is: (NIR - RED) / (NIR + RED)
class : RasterLayer
dimensions : 300, 300, 90000 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 390045, 399045, 4482105, 4491105 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : layer
values : -0.3114754, 0.5664336 (min, max)
# uncomment and run in 3D
plot3D(dem, drape = ndvi, zfac = 1.5)
### remove values below zero
ndvi[ndvi <= 0] <- NA
plot(ndvi)## plot again plot3D(dem, drape=ndvi)
## different way to plot in 2D
library(rasterVis)
levelplot(ndvi, col.regions = bpy.colors(100))Now you can see the areas of the highest reflectance and thus most healthy vegetation in November
—- skip to task 9 if short on time
Task 8: Create new information from satellite imagery: SAVI
SAVI stands for Soils Adjusted Vegetation Index, and this is another calibrated view of the ground.
### another index SAVI (soil adjusted vegetation index)
ndvi <- (n4 - n3)/(n4 + n3)
savi <- (n4 - n3)/((n4 + n3) * 0.25) # with L=1 -> similar to NDVI
### let´s compare visually
par(mfrow = c(1, 2))
plot(savi, main = "SAVI")
plot(ndvi, main = "NDVI")Task 9: Unsupervised Classification with k-means
We would like to isolate and better see the clusters of growth within our image. We will run kmeans function on a composite image in order to cluster data based on similarity or similar groups!
# first, let's select an image and make it into a brick including ndvi
data(nov2)
data(nov1)
n2 <- raster(nov2)
n1 <- raster(nov1)
ndviclass : RasterLayer
dimensions : 300, 300, 90000 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 390045, 399045, 4482105, 4491105 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : layer
values : -0.3114754, 0.5664336 (min, max)
integer(0)
# create a new composite brick out of the available data
myNewBrick <- brick(n4, n3, n2, n1, ndvi)
splom(myNewBrick)
Next, run the kmeans classification. Beware that the kmeans function
does not tolerate NA/INF/NaN and similar values. Our new brick should
not have any but in future classification remember that you need to get
around them, either by exclusion or substitution via mean values.
# Run kmeans classification on the values in your new brick Read on
# Thresholding here:
# https://rspatial.org/raster/rs/3-basicmath.html#vegetation-indices
ICE_df <- as.data.frame(myNewBrick)
set.seed(99)
cluster_ICE <- kmeans(ICE_df, 10) ### kmeans, with 4 clusters
str(cluster_ICE)List of 9
$ cluster : int [1:90000] 4 4 4 6 6 7 6 6 4 1 ...
$ centers : num [1:10, 1:5] 77.8 33.1 93.3 66.8 44.9 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:10] "1" "2" "3" "4" ...
.. ..$ : chr [1:5] "X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band4.geo.asc" "X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band3.geo.asc" "X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band2.geo.asc" "X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band1.geo.asc" ...
$ totss : num 20611557
$ withinss : num [1:10] 169596 158034 131141 241482 231449 ...
$ tot.withinss: num 1874040
$ betweenss : num 18737517
$ size : int [1:10] 4612 10833 1783 6804 16124 8862 16155 3728 4620 16479
$ iter : int 10
$ ifault : int 0
- attr(*, "class")= chr "kmeans"
# cluster_ICE <- cluster::clara(ICE_df, 4) ### another option, clara, with 4
# clusters
# convert cluster information into a raster for plotting
clusters <- raster(myNewBrick) ## create an empty raster with same extent than ICE
clusters <- setValues(clusters, cluster_ICE$cluster) # convert cluster values into raster
clustersclass : RasterLayer
dimensions : 300, 300, 90000 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 390045, 399045, 4482105, 4491105 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : layer
values : 1, 10 (min, max)
# uncomment to plot the clusters in 3D over the DEM
plot3D(dem, drape = clusters, col = c("white", "green", "blue", "yellow"))
# calculate the average spectral signature of 1-4 bands of growth
ICE_mean <- zonal(myNewBrick, clusters, fun = "mean")
ICE_mean # see the values for ndvi (layer) being most distinct zone X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band4.geo.asc
[1,] 1 77.78751
[2,] 2 33.07745
[3,] 3 93.29893
[4,] 4 66.75720
[5,] 5 44.89792
[6,] 6 57.19567
[7,] 7 50.56335
[8,] 8 56.26502
[9,] 9 47.30498
[10,] 10 39.66339
X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band3.geo.asc
[1,] 41.72029
[2,] 31.33093
[3,] 40.36960
[4,] 43.96326
[5,] 37.63799
[6,] 42.74362
[7,] 40.45453
[8,] 50.90504
[9,] 45.32706
[10,] 34.33989
X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band2.geo.asc
[1,] 45.40590
[2,] 34.81261
[3,] 46.12900
[4,] 45.14065
[5,] 38.55848
[6,] 42.28154
[7,] 39.87242
[8,] 48.89753
[9,] 44.43831
[10,] 36.50592
X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band1.geo.asc layer
[1,] 58.13964 0.30225062
[2,] 52.30970 0.02607170
[3,] 58.31408 0.39513993
[4,] 58.60553 0.20678500
[5,] 54.86039 0.08794067
[6,] 56.87599 0.14495138
[7,] 55.61634 0.11117963
[8,] 62.34576 0.04988818
[9,] 59.38939 0.02106308
[10,] 53.31756 0.07190780
Note that you have aggregated the final cluster raster by using the
focal() function using the mean function on
the four clusters identified by kmeans() as similar.
Task 10: What is the trend in de-/afforestation? - Individual tree crown segmentation
The ITC (Individual Tree Crowns) delineation approach finds
local maxima within imagery that contains subtle color
differences, such as the canopy image provided. The
itcIMG() function designates these maxima as tree tops,
then uses a decision tree method to grow individual crowns around the
local maxima.
The image we use is based on LiDAR (Light Detection and Ranging) in
xyz format. To get the segmentation going you will need either
rgeos with itcSegment version 07 or the
terra and current itcSegment packages
# install.packages('terra') install.packages('itcSegment')
library(terra)
library(itcSegment)
imgData <- rast("../data/imgData.tif")
plot(imgData)# Use the itcIMG() function to detect and grow the individual crowns
se <- itcIMG(imgData, epsg = 32632)
summary(se) ID X Y CA_m2
Min. : 1.00 Min. :676744 Min. :5091658 Min. : 6.48
1st Qu.: 30.25 1st Qu.:676754 1st Qu.:5091661 1st Qu.:34.12
Median : 59.50 Median :676776 Median :5091663 Median :43.74
Mean : 59.50 Mean :676782 Mean :5091663 Mean :45.66
3rd Qu.: 88.75 3rd Qu.:676802 3rd Qu.:5091666 3rd Qu.:56.17
Max. :118.00 Max. :676832 Max. :5091669 Max. :95.58
class : SpatVector
geometry : polygons
dimensions : 118, 4 (geometries, attributes)
extent : 676744.4, 676837.1, 5091659, 5091747 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632)
names : ID X Y CA_m2
type : <int> <num> <num> <num>
values : 1 6.768e+05 5.092e+06 59.94
2 6.768e+05 5.092e+06 82.62
3 6.768e+05 5.092e+06 54.27
### Let´s overlay the image and the product of segmentation (run both lines)
plot(imgData)
plot(se, axes = F, add = T)Task 11: Visualise the segmentation result in Leaflet
You can probably do all of this yourself, but here is a hint about projecting the SpatialPolygonsDataFrame, just in case:
class : SpatVector
geometry : polygons
dimensions : 118, 4 (geometries, attributes)
extent : 676744.4, 676837.1, 5091659, 5091747 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632)
names : ID X Y CA_m2
type : <int> <num> <num> <num>
values : 1 6.768e+05 5.092e+06 59.94
2 6.768e+05 5.092e+06 82.62
3 6.768e+05 5.092e+06 54.27
[1] "PROJCRS[\"WGS 84 / UTM zone 32N\",\n BASEGEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 32N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",9,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n AREA[\"Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],\n BBOX[0,6,84,12]],\n ID[\"EPSG\",32632]]"
# Project the SpatialPolygon using the spTransform() function se4326 <-
# spTransform(se,CRS( '+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# '))
# Project the SpatialPolygon using sf library
library(sf)
se4326 <- st_transform(st_as_sf(se), 4326)
# ...then we can combine them with leaflet
library(leaflet)
leaflet() %>%
addTiles() %>%
addRasterImage(imgData) %>%
addPolygons(data = se4326, weight = 1, color = "black") # TadaaCan you figure out where this forested landscape is from?
Task 12: OPTIONAL: a more demanding segmentation example with a bigger raster!
And because we have it in the data/ folder, try with this larger example.
r <- rast("../data/myDem_subset.tif")
r
plot(r)
r.se <- itcIMG(r, epsg = 25829, ischm = T) # can be a bit slow (2-3mins)
summary(r.se)
plot(r)
plot(r.se, axes = F, add = TRUE)
# Adjust 'th' argument for excessive capture of small growth
r.se5 <- itcIMG(r, epsg = 25829, th = 5, ischm = T) # th - how low should algorithm be looking for canopy
plot(r)
plot(r.se5, axes = T, add = TRUE)
# Write the result to shapefile
td
st_write(st_as_sf(r.se), "../data/itcTrees_subset", driver = "ESRI Shapefile")
# want to see it in Leaflet?
library(sf)
rse <- st_read("../data/itcTrees_subset/itcTrees_subset.shp")
rse4326 <- st_transform(rse, crs = 4326)
# Control question: where is this landscape from?
leaflet() %>%
addTiles() %>%
addProviderTiles("Esri.WorldPhysical") %>%
# addProviderTiles('Esri.WorldImagery') %>%
addRasterImage(r) %>%
addPolygons(data = rse4326, weight = 1, color = "black") # Neat :)The End
Similar approaches can be used when mapping socio-cultural phenomena in satellite imagery, from mosaicing images of night lights as proxies of economic performance, or detecting phenomena in the landscape such as burial mounds, growing urban sprawl, or tracing the outlines of scanned line drawings. (In the latter two you may need to base the classification on reflectance or edge detection rather than elevation.)
Bibliography
https://geoscripting-wur.github.io/AdvancedRasterAnalysis/ http://rspatial.org/spatial/rst/8-rastermanip.html http://neondataskills.org/R/Image-Raster-Data-In-R/ https://geoscripting-wur.github.io/IntroToRaster/ http://wiki.landscapetoolbox.org/doku.php/remote_sensing_methods:home https://rpubs.com/alobo/vectorOnraster